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A new near-to-far-field transformation algorithm for three-dimensional 
finite-different time-domain is presented in this article. This new approach is 
based directly on the polarization current of the scatterer, not the scattered 
near fields. It therefore eliminates the numerical errors originating from the 
spatial offset of the E and H fields, inherent in the standard near-to-far-field 
transformation. The proposed method is validated via direct comparisons 
with the analytical Lorentz-Mie solutions of plane waves scattered by large 
dielectric and metallic spheres with strong forward-scattering lobes. © 2009 
Optical Society of America 

The grid-based finite-difference time-domain (FDTD) method is one of the most popular 
Maxwell solvers, which has been proven to be efficient, stable and easy-to-implement [1]. Due 
to the computational resource limitation, a FDTD simulation truncates the open boundary 
to a spatial domain adjacent to the scatterer. The near-to-far-field (NTFF) transformation, 
therefore, is routinely employed to obtain the far-zone information such as antenna scattering 
patterns and nanocavity radiation patterns [1-8]. The standard NTFF (S-NTFF) transfor- 
mation, introduced in the early 1980s, is based on the vector Kirchhoff integral relation [1]. 
The scattered fields in the far zone are calculated through an integration of the near-zone 
fields over a virtual closed (Huygens) surface completely enclosing the scatterer [9]. To ac- 
complish this we need compute, via FDTD and discrete Fourier transformation, the scattered 
E and H fields tangential to the fictitious surface. However, the spatial and temporal offset 
between E and H field, a character of the Yee update scheme [1] , may result in unacceptable 
numerical errors [10-12]. For instance, Ref. [13] demonstrated that the accuracy of the S- 
NTFF is unacceptable when calculating the backscattering from strongly forward-scattering 
objects, and the relative error may be as high as two orders of magnitude at short wave- 
length. To improve its numerical accuracy, at least two different modifications have been 
proposed, including discarding information in the forward-scattering region [13] and using 
geometric mean in place of arithmetic mean [14]. 
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In classical electrodynamics, the optical response of medium made up of a large number 
of atoms or molecules originates from the perturbed motions of the charges bound in each 
molecule. The molecule charge density is distorted by the external electromagnetic fields, 
and further produces an electric polarization P in the medium [9]. Consequently in this 
article we proposal using the polarization current, rather than the scattered near fields, to 
derive the far-zone information. The proposed approach is validated via comparison with 
rigorous analytical solutions of plane wave scattered by large dielectric and metallic spheres. 
A detailed discussion regarding its advantages and disadvantages is also presented. 

In the absence of sources, the Maxwell equations are read as 

<9B 

V-B = 0, VxE=- — , 

at 
<9D 

V-D = 0, VxH = — . (1) 

Assuming the following constitutive relations D = eoE + P as well as B = //oH (the scatterer 
is therefore nonmagnetic), we arrive an inhomogeneous Helmholtz wave equation for the 
vector potential A in the Lorentz gauge [9] 

V 2 A-^ = -/i J(t), (2) 

where J is the polarization current defined as J = dP/dt [15]. With a time dependence e~ luJt 
understood, Eq. (2) becomes 

V 2 A + A; 2 A = -fi 3(u), (3) 

with k = uj/c being the vacuum wave number. A formal solution of the above equation, with 
the help of free-space Green function, can be written as [9] 
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The left-hand side represents the scattered wave with A^ ) representing the incident wave. To 
obtain the field in the radiation zone, it is sufficient to approximate the numerator with [9] 

|r — r'| ~ r — n ■ r', (5) 

while in the denominator |r — r'| ~ r. Here n is a unit vector in the direction of r. The 
far- zone vector potential of the scattered field is therefore expressed as 

lim A(r) « ^— ( J(rOe- lfcn r 'rfr' = ^^p. (6) 

It behaves as an outgoing spherical wave while depends on the polar and azimuth angles 
(9,(j>) in spherical polar coordinate. The scattered fields in the radiation zone, keeping only 
the leading order, are further given as 

B = !^ e ifcr n x p, E = cB x n. (7) 
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Fig. 1. Forward and backward scattering magnitudes of l-«m and 3-/im diameter dielectric 
spheres with e r = 1.21. Both numerical FDTD solution and rigorous Lorentz-Mie solution 
are presented. The cell size employed in FDTD simulation is 25 nm. The corresponding 
S-NTFF results can be found from Fig.l of Ref. [13]. 



Evidently, both fields are transverse to the radius vector n and fall off as r 1 . 
The time-averaged power scattered per unit solid angle is 



dP 

dn 



-Re 



r 2 n ■ E x H" 
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(n x p) x n| 



with rjo = y /in/en being the intrinsic impedance of free space. To obtain the scattering cross 
section, the scattered power needs be normalized to the incident power [9] . It is worth noting 
that the above general expression is quite similar to that of an oscillating electric or magnetic 
dipole, except that the p(n) here is angular-dependent [9]. 

The above derivation immediately implies two features of our approach: (1) It depends 
only on the polarization current J but not the scattered electromagnetic fields E and H, 
the error induced by the spatial offset between them in the S-NTFF is therefore completely 
avoided; (2) Unlike the S-NTFF, there is no requirement of scattered fields. Consequently, 
our approach is easier to implement for strongly convergent light sources such as Gaussian 
beam [16]. 

It should be emphasized that we can very easily extract the polarization current from 
FDTD, because an auxiliary differential equation of the current is frequently adopted by 
standard FDTD to simulate dispersive, nonlinear and even quantum- mechanical media [1]. 
For the special case of a dielectric scatterer, the standard FDTD itself is enough to provide 
the information of the current, since it connects to the electric field as J (a;) = — za>eoxE(a;) 
with x being the susceptibility. 

To validate our approach as well as the corresponding numerical algorithm, we begin by 
considering the backscattering from two strongly forward-scattering dielectric spheres with 
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a relative permittivity e r = 1.21 = 0.21) which were considered in Ref. [13,14] with the 
use of the (modified) S-NTFF. Assuming the incident field propagates along the positive x 
direction, the backscattering (corresponding to n = (—1,0,0)) power is simply given as 



fc 2 ?7o 
32tt 2 



b,(n)| 2 + |p 2 (n)| 2 , (9) 



that is, we do not need consider the x-component current. In our simulations, the FDTD 
lattice is terminated by a uniaxial perfectly matched layer (UMPL) in all directions [1], 
and the whole computational space is divided into uniform cubic cells with size of 25 ran. 
Furthermore, an impulsive wideband plane wave is excited to obtain the spectra in the 
wavelength domain 500-1000 ran. 

Fig.l presents the FDTD calculated far-field spectra for two spheres with diameters 
d = Iran and d = 3/xm. Because their sizes are on the order of incident wavelength, the 
appearance of higher-order multipoles distort the symmetrical electric-dipole pattern and 
result in strong forward scattering. The biggest contrast between forward and backward 
scattering appears for the 1/im-diameter sphere. At 500-ran wavelength, its scattering in the 
forward direction is almost three orders of magnitude stronger than its backward counter- 
part. As a reference, the analytical Lorentz-Mie solutions are also shown. Clearly, excellent 
agreement is achieved even when the incident wavelength is coarsely resolved on the FDTD 
grid. Under the same conditions, the relative error of the S-NTFF result is as high as two 
orders of magnitude [13], and modifications such as the geometric mean must be introduced 
to improve its numerical accuracy [14]. 

We further consider the scattering of a 1/im-diameter gold sphere with Drude-type per- 
mittivity approximated as 



e(o;) = 1.0- I (10) 

in the wavelength region 700-1200 ran. The bulk plasma frequency u p is taken as 1.367 x 
10 16 s _1 , and the phenomenological collision frequency 7 = 6.478 x 10 13 s _1 [17]. Because of 
the negative permittivity, surface plasmon polaritons (SPP) are excited and the resulting 
evanescent waves concentrate around the metallic surface. To simulate these nonpropagating 
fields, sufficiently finer resolution is needed in FDTD to achieve acceptable accuracy. Fig. 2 
shows the comparison of the analytical Lorentz-Mie solution with our FDTD numerical 
results (with cell size of 10 ran). Clearly, they are quite close, except tiny blue shifts of SPP 
resonances as well as a relative difference of roughly 5% around the 770-nm wavelength. It 
should be emphasized that these disagreements are induced by the inherent staircasing of 
the standard FDTD on a Cartesian grid and not due to our NTFF. In the simulation the 
diameter of the metallic sphere is resolved by 100 cells, and better agreement can be achieved 
by decreasing the cell size. 
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Fig. 2. Far-field scattering spectra of 1 urn-diameter gold sphere. Both numerical FDTD 
solution and rigorous Lorentz-Mie solution are presented. The cell size employed in FDTD 
simulation is 10 nm. 

The key differences between our current-based NTFF with the S-NTFF can be summa- 
rized as follow: (1) Origins of numerical errors are different. They both suffer from staircasing. 
However, the S-NTFF has two more sources including the scattered-field requirement (we 
thus need total field/scattered field or pure scattered field scheme [1]) and the spatial offset 
of the E and H fields over the integral surface [13,14]; (2) A virtual integral surface enclos- 
ing the whole scatterers is not needed in our approach. It is therefore easily and efficiently 
implemented within the framework of standard FDTD codes; (3) The two methods have 
different computational burdens. Consider a cubic scatterer resolved by N 3 cells, the total 
number of discrete Fourier transformations to be performed is 3 x N 3 in volume integration 
(our method) in comparison with 12 x N 2 in surface integration (S-NTFF). On the contrary, 
a multiple scatterer cluster with small volume but large surface may require more calcu- 
lational resources in the S-NTFF; (4) The current-based NTFF has a close relation with 
the analytical multipole expansion, it is therefore convenient to identify the contributions 
from different multipolar sources and further provide insights into complex electromagnetic 
processes including radiation, diffraction and scattering [5,9]. 

It should be pointed out that the spatial offset between E and H field can be avoided by 
using the scalar Kirchhoff integral formula 



with R = r — r' and ip being any component of E or H [9]. However, the integral surface 
S generally is infinite for an open-boundary problem. To achieve acceptable accuracy in 
numerical simulation, we need to truncate it to a finite surface with sufficient area. Consider 
an electric dipole with size a. The radius of the integral surface should be roughly 10a to 
contain the field with amplitude bigger than 10% of the dipole magnitude. It thus increases 
the computational burden. 

To conclude, a polarization-current-based near-to-far-field transformation for the three- 




(11) 
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dimensional finite-difference time-domain algorithm is developed. Because it completely elim- 
inates the spatial-offset error inherent in the surface integration of standard near-to-far-field 
scheme, the new method can achieve improved accuracy even with coarse grids. Its valid- 
ity and efficiency is further demonstrated via direct comparisons with analytical theory by 
calculating the backscattering from strongly forward-scattering objects including dielectric 
and metallic spheres. We also compare it to the standard near-to-far-field transformation in 
terms of calculational cost and accuracy. The new method may have potential applications 
in biophotonics and nanophotonics. 
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